{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lab 11 - Logistic Regression Continued\n", "\n", "The Akimel O'odham people, who were also known as the Pima Indians since European colonization of the US, currently have a high prevalence of diabetes. The Pima Indian Diabetes dataset contains different possible diabetes indicators and whether the person has diabetes is on [Kaggle](ttps://www.kaggle.com/uciml/pima-indians-diabetes-database) or available [here](http://comet.lehman.cuny.edu/owen/teaching/mat328/diabetes.csv).\n", "\n", "Load the dataset into a dataframe called `diabetes`." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import numpy as np\n", "from sklearn.model_selection import train_test_split\n", "from sklearn.linear_model import LogisticRegression\n", "import statsmodels.formula.api as smf\n", "from scipy.special import expit\n", "from scipy.stats import logistic\n", "\n", "%matplotlib inline\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
PregnanciesGlucoseBloodPressureSkinThicknessInsulinBMIDiabetesPedigreeFunctionAgeOutcome
061487235033.60.627501
11856629026.60.351310
28183640023.30.672321
318966239428.10.167210
40137403516843.12.288331
\n", "
" ], "text/plain": [ " Pregnancies Glucose BloodPressure SkinThickness Insulin BMI \\\n", "0 6 148 72 35 0 33.6 \n", "1 1 85 66 29 0 26.6 \n", "2 8 183 64 0 0 23.3 \n", "3 1 89 66 23 94 28.1 \n", "4 0 137 40 35 168 43.1 \n", "\n", " DiabetesPedigreeFunction Age Outcome \n", "0 0.627 50 1 \n", "1 0.351 31 0 \n", "2 0.672 32 1 \n", "3 0.167 21 0 \n", "4 2.288 33 1 " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "diabetes = pd.read_csv(\"diabetes.csv\")\n", "diabetes.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot a scatter plot of glucose vs. diabetes outcome." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmYHHd95/H3t6qvuceaGVmnpdGFD2yCmRixJg4LxNiG4NxrNsQhJniTBUI2hKx5yBIedrMJYQMJrLOJeWKD/YQQ2FxKAmuIIXE4bDzC+JItocOydVkzOubqmemjvvtH1ZRb4zll9bQsfV7PM4+6fv3rqq+qq+rTdXS1uTsiIiIAQaMLEBGRs4dCQUREUgoFERFJKRRERCSlUBARkZRCQUREUgoFERFJKRRERCSlUBARkVSm0QUsVnd3t69fv77RZYiIvKRs37590N175uv3kguF9evX09/f3+gyREReUsxs/0L66fCRiIikFAoiIpJSKIiISEqhICIiKYWCiIikFAoiIpJSKIiISKpu31MwszuBtwBH3f3lMzxvwB8DNwBF4B3u/r161TM4OsGhkxOs6izQ3Vqo12SkQYqlCsPjFdqb4kV66nFzLnNKn+eGJ8CNCzvy6XO1ry2WKuly0pzLMDxeIRNCpQrtTRmOjU6y++gom5a30pQLT1mmdh0Z5vGDQ2xc3srYZJn9x8e5dGU71QiOjkxwQXOWE2Nllrfnac1nePiZk5g5nU1ZDpyc4FXrLmBZS47dR0cJDPYNFuluybFlRRv9+4/z9ECRqzd30dNWYHSiSmdzhl1HRjhWLLG+q5kTxQotuYDl7QWePjZGczYklwn4yuOHGZ+s8sZLl7Ouu41VnQX2HB3lO3uO0dGU4Zlj4yxvz3HlumUcGRqnOZdhXVdLOu6LV7QxMDzB4GiZS1a1MTxeZqxUobe7laNDk4yVK/R2tzA0XgaMFe15JioR46UK258+jgMXLWti19ExLmjKcdnqdk6MlVjeXuCirmaeOVbk6PAEF7TkKJaqbFreSldrnu8/c4L9x4qs6iwQOWxa3sqhk+Pcv2uAl69q54c3dFGpQrla5bt7Bzl0coKXrWhjvBLRnA25eGU7IxMVRieq5LPGYwdOMl6uclVvF9kwSNsPnxxneXuBcjXiwb3HWd1Z4KoNXelrV3bmyYYhmRCeOTbG0ZESTdmA7+0/yWSlyr+/eDlthWw6vuNjpRf8Hzb0tDA8XqZYqnLJynbam3Lp8jm1/JWrVU4UyxQyAfuPjVEsVVjR0cSJYpnlbTku6mpJl8Pa5fpMs3r9RrOZXQOMAnfPEgo3AO8lDoVXA3/s7q+eb7x9fX2+2C+v/dOjh/jk13ZRiZxMYLz/2i1cf/mqRY1Dzl5PHhrmngf3U6lGjEyUMYzWQoZMGHDz1nVcvLKdJw8N8+mv/4CnjowA8LIVrbzvDVtwJ33toRPj7D9eJAyMauSsW9ZMcz5k78AYG3taOToywc4jI1hgVCpV2gpZWgtZMoFx0QUF/m3PcaqRE82wShkw1RwYM/YByARxx8ocq6UBuRAmq4ufV7kAwsAYn2sCp8kAs/jf0KAUzd0/NMhnAibLEQ5EQMYglwloyYUMjpXTeZYxqPrz8xAgH8DLVrTx+KERZppUaJANjUrVXzA/s4FhBqWqn/LeTAksfi8yQYCZsbG7hYMnipycrFCdYWIBkAmMcuRkw/h1TdmA4+MVpm9iQ4OXr27nstWdvHZjF9/cc4x9A6M8cWgYPGK0FJ2yfFjyd0FLlr51y+hqy6fL9WKY2XZ375uvX90OH7n7/cDxObrcSBwY7u4PAJ1mtvJM1zE4OsEnv7aLQjbkwvYChWzIH351F4OjE2d6UtIAxVKFex7cT3M2pLs1z77BInsHR+luzdOcDbn7gf0Mjk5w17f3sXdwlI6mLJ3NWfYNFrnj/j3c9e19NGdDOpqyPH5oiNHJCstacoxOVnjs4EmeOT5OIRuyb3CUxw8OUY2cllxAOYJjY2U6m7IY8PVdx+It/Szb2trm2QIBoBLNHQhT4zqdQIB4Q12PQIC4rsjjjfd8gQBxv2I5wo10o15xKFUjBsbKp/StzDBrJyN4dJZAmBr/ZOWFgQBQjpxS1dO6p4scylWoViPK1So7jgxzfLxCNMvEIqAUOQ5UI5goRxwrVmYcedVhx+ERJstVPvG1XVSrEU8dGSHEXxAIaX0Gx8fK7Ds2Si407n5gP8VSZZb/+YvTyHMKq4Fna4YPJG0vYGa3mlm/mfUPDAwsaiKHTk5QiZyWfLy71ZLPUImcQycVCueC4fEKlWpESz7DZDkiMAjMmKzEbZVqxKGTExRLVQIzcpmAbBgQGAxNVCiWqrTkM4xMxCtYGBiT5YgwsHjjW4mfH08+zYaBUak+/6l4vFylmnwUDEObLRNkDtM/SU8N25kY94sdgUFgQRoG89VkxPVP9bNZXhC5MzJRphI54+UqkTuZTPiCeTF93OOliDAIqFQjhsfPvVBYMHe/w9373L2vp2fe+zmdYlVngUxgjE3GM3BsskImMFZ16rzCuaC9KT5MNDZZIZ8NiDxe4fKZuC0TBsn5gZDInVIlolyNP411FDI050LGJiu0FeIPDdXIyWcDqlF8WCGfiZ9vygZY8nwmjDc27tCUDQmTNb+aHIqQxZm+4ZwaPhMB+6LfD4fII4IgHZyvO2bP95ttIx+Y0ZYcemzKhgQWH5KcLUSmxt2UC6hGEZkwSM+fnWmNDIWDwNqa4TVJ2xnV3Vrg/dduYaJc5bnhCSbKVd5/7RadbD5HNOcy3Lx1HcVylcHRSXq7m9nQ3crg6CTFcpWbt66ju7XALVf3sqG7laHxMieLZXq7m7n1mo3ccnUvxXKVofEyL1/VQWs+w/GxEq35DJev7uSiZU1MlKv0drdy+eoOwsAYK0VkA+hqyXJyPD7u/YYtXfGB6FlW6trmYI4VPxPEx8/nYkA+XOSMSuQCaJpvAqfJiP9vocXTmU9o0JwNMH9+Q5QxyIUBPS3ZU/pmZpi1+QCuWNU260YsY5DP2IzzMxsYudDSuqcLDLIhhGFANgy5dEU7y5oyaTi8oD+QCyw+nxJAIRvQ1ZyZceShwaUr28hnQ95/7RbCMODiFW1UMVpzwQuWDwNwWNaSpberlVLVuXnrurqdbK7biWYAM1sP/OMsJ5rfDLyH5080f8rdr5pvnKdzohl09dG5Tlcf6eojXX00t4WeaK7n1Ud/CbwO6AaeA34HyAK4+58ml6T+b+A64ktSf8nd593an24oiIiczxYaCnW72NXd3zbP8w68u17TFxGRxXtJnGgWEZGloVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkVddQMLPrzGynme02s9tmeP4iM/uGmT1sZo+a2Q31rEdEROZWt1AwsxC4HbgeuBR4m5ldOq3bbwNfdPdXAjcBf1KvekREZH713FO4Ctjt7nvdvQR8AbhxWh8H2pPHHcChOtYjIiLzyNRx3KuBZ2uGDwCvntbnI8BXzey9QAvwxjrWIyIi82j0iea3AZ919zXADcA9ZvaCmszsVjPrN7P+gYGBJS9SROR8Uc9QOAisrRlek7TVeifwRQB3/w5QALqnj8jd73D3Pnfv6+npqVO5IiJSz1B4CNhsZr1mliM+kbxtWp9ngDcAmNklxKGgXQERkQapWyi4ewV4D3Av8CTxVUZPmNlHzeytSbf3A+8ys0eAvwTe4e5er5pERGRu9TzRjLt/GfjytLYP1zzeAVxdzxpERGThGn2iWUREziIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSdU1FMzsOjPbaWa7zey2Wfr8nJntMLMnzOzz9axHRETmlqnXiM0sBG4Hfgw4ADxkZtvcfUdNn83AB4Gr3f2EmS2vVz0iIjK/eu4pXAXsdve97l4CvgDcOK3Pu4Db3f0EgLsfrWM9IiIyjwWFgpk1m9l/M7PPJMObzewt87xsNfBszfCBpK3WFmCLmX3LzB4ws+sWWriIiJx5C91TuAuYBF6TDB8E/scZmH4G2Ay8Dngb8Bkz65zeycxuNbN+M+sfGBg4A5MVEZGZLDQUNrr7HwBlAHcvAjbPaw4Ca2uG1yRttQ4A29y97O77gF3EIXEKd7/D3fvcva+np2eBJYuIyGItNBRKZtYEOICZbSTec5jLQ8BmM+s1sxxwE7BtWp+/I95LwMy6iQ8n7V1gTSIicoYt9Oqj3wH+H7DWzP4CuBp4x1wvcPeKmb0HuBcIgTvd/Qkz+yjQ7+7bkueuNbMdQBX4gLsfO73/ioiIvFjm7gvraNYFbCU+bPSAuw/Ws7DZ9PX1eX9/fyMmLSLykmVm2929b75+i7kkdTXxJ/4ccI2Z/dTpFiciImenBR0+MrM7gSuAJ4AoaXbgb+pUl4iINMBCzylsdfdL61qJiIg03EIPH33HzBQKIiLnuIXuKdxNHAxHiC9FNcDd/Yq6VSYiIktuoaHw58AvAI/x/DkFERE5xyw0FAaS7xWIiMg5bKGh8HDyWwf/QM03md1dVx+JiJxDFhoKTcRhcG1Nmy5JFRE5xywoFNz9l+pdiIiINN5Cf09hjZn9rZkdTf7+2szW1Ls4ERFZWov5PYVtwKrk7x+SNhEROYcsNBR63P0ud68kf58F9MMGIiLnmIWGwjEze7uZhcnf2wHd4lpE5Byz0FC4Bfg54AhwGPgZQCefRUTOMQu9+mg/8NY61yIiIg220KuPPmdmnTXDFyS30xYRkXPIQg8fXeHuJ6cG3P0E8Mr6lCQiIo2y0FAIzOyCqQEzW8bCvw0tIiIvEQvdsP8h8a2zv5QM/yzwP+tTkoiINMpCTzTfbWb9wOuTpp9y9x31K0tERBphob/RfI+7/wKwY4Y2ERE5Ryz0nMJltQNmFgKvOvPliIhII80ZCmb2QTMbAa4ws2EzG0mGjwJ/vyQViojIkpkzFNz999y9Dfi4u7e7e1vy1+XuH1yiGkVEZIks9Oqjr5jZNdMb3f3+M1yPiIg00EJD4QM1jwvAVcB2nr8aSUREzgELvST1x2uHzWwt8Ed1qUhERBpmoVcfTXcAuORMFiIiIo230O8pfBrwZDAgvu/R9+pVlIiINMZC9xR2ALuSvweA33L3t8/3IjO7zsx2mtluM7ttjn4/bWZuZn0LrEdEROpgzj0FM8sQ3+PoFuCZpPki4E4z+667l+d4bQjcDvwY8eGmh8xs2/TbY5hZG/A+4MHT/l+IiMgZMd+ewseBZUCvu1/p7lcCG4BO4H/N89qrgN3uvtfdS8AXgBtn6PffgY8BE4uqXEREzrj5QuEtwLvcfWSqwd2HgV8FbpjntauBZ2uGDyRtKTO7Eljr7v8014jM7FYz6zez/oGBgXkmKyIip2u+UHB39xkaqzx/4vm0mFkAfAJ4/3x93f0Od+9z976enp4XM1kREZnDfKGww8xunt5oZm8HnprntQeBtTXDa5K2KW3Ay4F/MbOnga3ANp1sFhFpnPkuSX038DdmdgvxN5gB+oAm4Cfnee1DwGYz6yUOg5uA/zj1pLsPAd1Tw2b2L8Bvunv/Yv4DIiJy5swZCu5+EHi1mb2e52+f/WV3v2++Ebt7xczeA9wLhMCd7v6EmX0U6Hf3bS+ydhEROcNshlMGZ7W+vj7v79fOhIjIYpjZdnef9/D86d7mQkREzkEKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIKBRERSdU1FMzsOjPbaWa7zey2GZ7/DTPbYWaPmtl9ZraunvWIiMjc6hYKZhYCtwPXA5cCbzOzS6d1exjoc/crgP8L/EG96hERkfnVc0/hKmC3u+919xLwBeDG2g7u/g13LyaDDwBr6liPiIjMo56hsBp4tmb4QNI2m3cCX5npCTO71cz6zax/YGDgDJYoIiK1zooTzWb2dqAP+PhMz7v7He7e5+59PT09S1uciMh5JFPHcR8E1tYMr0naTmFmbwQ+BPyou0/WsR4REZlHPfcUHgI2m1mvmeWAm4BttR3M7JXAnwFvdfejdaxFREQWoG6h4O4V4D3AvcCTwBfd/Qkz+6iZvTXp9nGgFfiSmX3fzLbNMjoREVkC9Tx8hLt/GfjytLYP1zx+Yz2nLyIii3NWnGgWEZGzg0JBRERSCgUREUkpFEREJKVQEBGRlEJBRERSCgUREUkpFEREJKVQEBGRlEJBRERSCgUREUkpFEREJKVQEBGRlEJBRERSCgUREUkpFEREJKVQEBGRlEJBRERSCgUREUkpFEREJKVQEBGRlEJBRERSCgUREUkpFEREJKVQEBGRlEJBRERSCgUREUkpFEREJKVQEBGRVF1DwcyuM7OdZrbbzG6b4fm8mf1V8vyDZra+nvWIiMjcMvUasZmFwO3AjwEHgIfMbJu776jp9k7ghLtvMrObgI8B/6Ee9Tx7fIzdR0fZtLyVtcta6jGJl5RiqcLweIX2pgzNudNfDIqlCk8PjjI6WWVDTwvdrYUXPP/c8AS4cWFHPp1W7fQBnjo8zNHhSS5b3Z6+P7uODLN9/3HWdbXQ3Zrnu3uPUSxHvGJtB025DIVMwJGhSTCnoynLziPDAHQ2ZfnW7kEmKxFXb+ri4Mlxjo6UaM2FbN9/grZ8yJtfsYq9A0UOnCiysqPA44eGaMmG9Pa08s1dA3hgvOqiTsYmI1Z1Ftg9MEq5GrGqo0D/vhNks8ZrNnRzZGSCNR0FDpwo8pXHj5AN4XVblnNouEQ1cl6xpoPnhkqs72mmNZ/hi9/dz2TkXLWuk2dPTFCqOj+8/gIgoKctx4N7BxkYK3HJha0cG6+yqbuFFR0F7t81wOrOJn7ook7u2/EcQ5NlfnTzcgZHyhwrjrOms5ldz42w8oImVrQVuO+pozTnQl6zoYsnDw9RrkRcvLKDofEKqy8ocPBEEcdY3pbn4WdOsqK9wPVXrGR4vMzYZMTKzjxHhidozoas62ph15EhBscqtOYDdh4Zpbe7me7WPP37T9Db1cylqzvZeXgEN2dFe4FjoyW6WnOMlyOWtxUIA/jmrgEMuHpLDxe2N9HelKFYqvDw/pMUyxUuWdFOJXJGJ6qMl8s8eXiELRe2snVjN8VShUMnJyhkAk6Ol2nNZ1nf3UxzLh7Hc8lycGF7vPwNj1cYHi+xd3CM5W05Ll7ZMe9yPn2ZnGn9GByd4NDJCS5ozpINw/T5haxPs41/6nEmhEqVU8ZxptbTxTB3r8+IzV4DfMTd35QMfxDA3X+vps+9SZ/vmFkGOAL0+BxF9fX1eX9//6Jq+bN/3c3t39hD1Z3QjF97/UZ++ZpNp/G/Ojc8eWiYex7cT6UakQkDbt66jotXtp/WeH73yzt45NkhHKezKcuH3nwJ11++Kn3+01//AU8dGQHgZStaed8btuBOOv3RiQq7jw6z99g4OOQzAf/ljZs4PDTBZ7+9n8hhpoUhH0K5CmZQrc8iLGeI8cL3cHN3Mxd2NPHowZOMTFTT53OhUY38lPd0WVNASyHPeLnCyWKFTADZTMgVq9v5+a3r+IdHDrHzyCgAqzoLtBUyHDwxzo7DI0TuBGa8fFUbv//Tr5h1Oa9dJ0YmyhhGayFzyvrxT48e4pNf20WxXGV0osJlK9voXd7Gazd18c3dx+Zcn2rHPzpRwXHaCtn0cTVy9g6MsbGnla62PDdvXXfKevJi1tP0fTDb7u598/Wr5+Gj1cCzNcMHkrYZ+7h7BRgCus5kEc8eH+P2b+whGxidTVmygfGpr+/h2eNjZ3IyLxnFUoV7HtxPczZkZUcTzdmQux/YT7FUWfR4PvPNPTx2cIhCNqA1n2GsVOXj9+5kcHSCYqnCXd/ex97BUTqasnQ2Z9k3WOSO+/dw17f30ZwN6W7N89SRYXYPFMng5DMBlWrEJ/95N3d+az/uENjM05+sQoQC4aVgprdo92CR7fuPMzxRpfYtLlX9Be/p8fGIY6PjjExUiSKnXHVygfHYwWF+/ytPsmdgjM7mLG2FDI8dHOKpwyM8eWSEKHICg0wATxwe4U/+5QczLue160R3a559g0X2Do7S3ZpP149nj4/xya/tIpcJKFcjsoHx1HOjVCpVPvHVXeRCm3V9mj7+vYOj7Bss0tGUZe/gKHuOjnHgxDiFbMjBoXFyoXHnt/al68mLWU9Px0viRLOZ3Wpm/WbWPzAwsKjX7j46StWdQi4EoJALqbqz++hoPUo96w2PV6hUI1ry8a5oSz5DpRoxPL64hS3ePY9fk8sEhIERBkapEnHo5ATD4xWKpSqBGblMQDYMCAyGJuL2lnyGyUpEueq4Q5g8HwRGuRrhQPCSWDrldDhQjeKt/2zBX6saEW/kk2UiSF40UYrASZavuG2yEuE1fcMgAIfjY+UZl/PadWKyHMXLoRmTlShdP3YfHaUSOflMiHu8HXF3JsoRlcjjaTDz+nTK+CsRgRmBwchEhcAMx6lUnZZ8higZV7FUTdeT2cZbL/Vc7Q4Ca2uG1yRtM/ZJDh91AMemj8jd73D3Pnfv6+npWVQRm5a3EpoxUaoCMFGqEpqxaXnrosZzrmhvineJxybjhWtsskImDNJjm4sZz9RrSpUo3uWPnFwmYFVnITkGGhK5U6pElKsRkUNHIW4fm6yQzwRkQ4sPASXPR5GTDQMMiKIz/b+Xs4UBYbJhjxawtxcGcRBMLRNR8qJCLgAjWb7itnwmwGr6VqMIDJa1ZGdczmvXiXw2iJdDj/dcp9aPTctbyQTGZKWKWbwdMTMK2YBMYPE0mHl9OmX8mYDIncihrZAhcscwMqExNlkhSMbVnAvT9WS28dZLPUPhIWCzmfWaWQ64Cdg2rc824BeTxz8DfH2u8wmnY+2yFn7t9RspR87J8TLlyPm11288b082N+cy3Lx1HcVylcND4xTLVW7eum7RJ7Gacxlu/ZGNXL66g4lyxOhkhZZcyAfe9DK6Wws05zLccnUvG7pbGRovc7JYpre7mVuv2cgtV/dSLFcZHJ3k4hXtbOpppkL8ySwTxucUfvm16zCbfYORD+OFN1zAp0xprJneok3dzbxq3TLaC+Eph5dyob3gPV3WFNDV2kRbISQIjGxolCLn8tXtfPD6S9jY08LJYpmRiQqXr+7g4pVtXLKiLQ4Rh0oEl61s4z+/bvOMy3ntOjE4OklvdzMbulsZHJ1M14+1y1p4/7VbKFUismFAOXIuvrCVTCaM26s+6/o0ffwbulvp7W5maLzMhu5WNi5vYc0FTUyUq6zuaKJUdW65ujddT17Meno66naiGcDMbgD+CAiBO939d83so0C/u28zswJwD/BK4Dhwk7vvnWucp3OiGXT10XS6+khXH+nqo9nXialxnEtXHy30RHNdQ6EeTjcURETOZ2fD1UciIvISo1AQEZGUQkFERFIKBRERSSkUREQkpVAQEZGUQkFERFIvue8pmNkAsP80X94NDJ7Bcs6ks7U21bU4qmvxztbazrW61rn7vPcJesmFwothZv0L+fJGI5yttamuxVFdi3e21na+1qXDRyIiklIoiIhI6nwLhTsaXcAcztbaVNfiqK7FO1trOy/rOq/OKYiIyNzOtz0FERGZw3kTCmZ2nZntNLPdZnZbA+tYa2bfMLMdZvaEmb0vaf+ImR00s+8nfzc0oLanzeyxZPr9SdsyM/uamf0g+feCJa7pZTXz5PtmNmxmv96o+WVmd5rZUTN7vKZtxnlksU8ly9yjZnblEtf1cTN7Kpn235pZZ9K+3szGa+bdny5xXbO+d2b2wWR+7TSzN9Wrrjlq+6uaup42s+8n7Usyz+bYPizdMubu5/wf8Y/87AE2ADngEeDSBtWyErgyedwG7AIuBT4C/GaD59PTQPe0tj8Abkse3wZ8rMHv4xFgXaPmF3ANcCXw+HzzCLgB+Arxj49tBR5c4rquBTLJ44/V1LW+tl8D5teM712yHjwC5IHeZJ0Nl7K2ac//IfDhpZxnc2wflmwZO1/2FK4Cdrv7XncvAV8AbmxEIe5+2N2/lzweAZ4EVjeilgW6Efhc8vhzwE80sJY3AHvc/XS/vPiiufv9xL8SWGu2eXQjcLfHHgA6zWzlUtXl7l9196lfen+A+HfSl9Qs82s2NwJfcPdJd98H7CZed5e8NjMz4OeAv6zX9Gepabbtw5ItY+dLKKwGnq0ZPsBZsCE2s/XEP0X6YNL0nmQX8M6lPkyTcOCrZrbdzG5N2i5098PJ4yPAhQ2oa8pNnLqSNnp+TZltHp1Ny90txJ8op/Sa2cNm9q9m9iMNqGem9+5sml8/Ajzn7j+oaVvSeTZt+7Bky9j5EgpnHTNrBf4a+HV3Hwb+D7AR+CHgMPGu61J7rbtfCVwPvNvMrql90uP91YZcrmZmOeCtwJeSprNhfr1AI+fRbMzsQ0AF+Iuk6TBwkbu/EvgN4PNm1r6EJZ2V7900b+PUDyBLOs9m2D6k6r2MnS+hcBBYWzO8JmlrCDPLEr/hf+HufwPg7s+5e9XdI+Az1HG3eTbufjD59yjwt0kNz03tjib/Hl3quhLXA99z9+eSGhs+v2rMNo8avtyZ2TuAtwA/n2xMSA7PHEsebyc+dr9lqWqa471r+PwCMLMM8FPAX021LeU8m2n7wBIuY+dLKDwEbDaz3uQT503AtkYUkhyr/HPha2CrAAADXklEQVTgSXf/RE177XHAnwQen/7aOtfVYmZtU4+JT1I+TjyffjHp9ovA3y9lXTVO+eTW6Pk1zWzzaBtwc3KFyFZgqOYQQN2Z2XXAbwFvdfdiTXuPmYXJ4w3AZmDvEtY123u3DbjJzPJm1pvU9d2lqqvGG4Gn3P3AVMNSzbPZtg8s5TJW77PpZ8sf8Vn6XcQJ/6EG1vFa4l2/R4HvJ383APcAjyXt24CVS1zXBuIrPx4BnpiaR0AXcB/wA+CfgWUNmGctwDGgo6atIfOLOJgOA2Xi47fvnG0eEV8RcnuyzD0G9C1xXbuJjzdPLWd/mvT96eQ9/j7wPeDHl7iuWd874EPJ/NoJXL/U72XS/lngV6b1XZJ5Nsf2YcmWMX2jWUREUufL4SMREVkAhYKIiKQUCiIiklIoiIhISqEgIiIphYKc98zsQjP7vJntTW7x8R0z+0kze52Z/WOj6xNZSgoFOa8lXxb6O+B+d9/g7q8i/nLjkt88TuRsoFCQ893rgZK7p/fHd/f97v7p2k7JbwD8Zs3w48kNyzCzm5Obuz1iZvckbevN7OtJ+31mdlHS/rPJax8xs/uTttDi3z54KOn/n+r+vxaZRabRBYg02GXE31A9LWZ2GfDbwL9z90EzW5Y89Wngc+7+OTO7BfgU8e2OPwy8yd0PWvKjN8Tf8h1y9x82szzwLTP7qse3jxZZUtpTEKlhZrcnn+IfWuBLXg98yd0HAdx96v78rwE+nzy+h/j2BQDfAj5rZu8i/tEgiO8zdbPFv/L1IPEtDTa/uP+JyOnRnoKc754gvq8NAO7+bjPrBvqn9atw6oeowulMzN1/xcxeDbwZ2G5mryK+f8173f3e0xmnyJmkPQU5330dKJjZr9a0Nc/Q72nin24k+R3c3prX/6yZdSXPTR0++jbxCWuAnwf+LXl+o7s/6O4fBgaIb3t8L/CryS2TMbMtyZ1qRZac9hTkvObubmY/AXzSzH6LeEM9BvzXaV3/mvgQzxPEh3h2Ja9/wsx+F/hXM6sCDwPvAN4L3GVmH0jG+UvJeD5uZpuJ9w7uI74r7aPEvwH8veRqqAEa+7Onch7TXVJFRCSlw0ciIpJSKIiISEqhICIiKYWCiIikFAoiIpJSKIiISEqhICIiKYWCiIik/j/dA8l5xVWbXQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "diabetes.plot.scatter(x = \"Glucose\", y = \"Outcome\", alpha = 0.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit a logistic regression model to this data, using Glucose as the independent variable and Outcome as the dependent variable." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: 0.526510\n", " Iterations 6\n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Logit Regression Results
Dep. Variable: Outcome No. Observations: 768
Model: Logit Df Residuals: 766
Method: MLE Df Model: 1
Date: Thu, 10 Oct 2019 Pseudo R-squ.: 0.1860
Time: 17:40:58 Log-Likelihood: -404.36
converged: True LL-Null: -496.74
LLR p-value: 4.418e-42
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err z P>|z| [0.025 0.975]
Intercept -5.3501 0.421 -12.713 0.000 -6.175 -4.525
Glucose 0.0379 0.003 11.647 0.000 0.031 0.044
" ], "text/plain": [ "\n", "\"\"\"\n", " Logit Regression Results \n", "==============================================================================\n", "Dep. Variable: Outcome No. Observations: 768\n", "Model: Logit Df Residuals: 766\n", "Method: MLE Df Model: 1\n", "Date: Thu, 10 Oct 2019 Pseudo R-squ.: 0.1860\n", "Time: 17:40:58 Log-Likelihood: -404.36\n", "converged: True LL-Null: -496.74\n", " LLR p-value: 4.418e-42\n", "==============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept -5.3501 0.421 -12.713 0.000 -6.175 -4.525\n", "Glucose 0.0379 0.003 11.647 0.000 0.031 0.044\n", "==============================================================================\n", "\"\"\"" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "logit_model = smf.logit(\"Outcome ~ Glucose\", diabetes).fit()\n", "logit_model.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What is the equation of the logistic regression model?\n", "\n", "$$y = \\frac{1}{1 + e^{-(-5.3501 + 0.0379x)}}$$\n", "\n", "\n", "We can also plot the logistic regression model using Seaborn's `regplot()`. Use `regplot()` as if you were doing linear regression on the variables, but add in the parameter `logistic = True`." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8HPV9+P/Xe/bULVmyfBvbYGOMTTgcjoQAIUCAAAaagxzN2ULza/JNm5CUhISkJDQHaZukTVNoSnM0KS0Jh5OQQICAE4iJzWVsYxtjbGz5kCzJulZ7zbx/f8xqvZJW0krW6nw/H499aHd2dvazo9l573yO90dUFWOMMQbAGe8CGGOMmTgsKBhjjMmyoGCMMSbLgoIxxpgsCwrGGGOyLCgYY4zJsqBgjDEmy4KCMcaYLAsKxhhjsoLjXYDhqqur00WLFo13MYwxZlJ55plnDqvqzKHWm3RBYdGiRWzcuHG8i2GMMZOKiOwpZD2rPjLGGJNlQcEYY0yWBQVjjDFZFhSMMcZkWVAwxhiTVbTeRyJyF3AF0KiqK/M8L8C3gcuBGPBBVX22GGX5ziM7+P4fXqUr6VIWDvAX5y7m/120rBhvZcbB49sauWPdLva2xigPBxAROhJpFtSUcsN5S7hgeT2Pb2vk67/Zxq7DXQAsri3lpstOAsi+FlVaulIkXI9IwGFGWYjulEvKVcJB//dTc2eSnmmpQgHBUygLB6iMBNjXlhiPjz9sAkyUqbXylUWAgCOkvaFLOdhnGa3PKYAIqBa+PcncvAGeLw0HeN38as5ZMoM/7mphy/422uPpAbfnCIQCDkvqyvi7S5dzwfL64X2IYZBizbwmIucBncCPBggKlwMfxw8KZwHfVtWzhtru6tWrdThdUr/zyA6+/dhOHPF3rKf+7RMXnmCBYQp4fFsjt6zdQiggpF2PhiNxAOZVRwkGHFKu8vbT5/Hj9XtojaVwxH+dp1AacgiHAlSVhGjvTtLUmQIgIOBmvhaO+CeolJv/e+Iw8BffmMFEgw4pTykJCp3JoY+inmBZUxri9re/btiBQUSeUdXVQ61XtOojVV0HtAyyyhr8gKGquh6oFpE5o12O7//hVRyBoOPgiJP56y83k98d63YRCgil4SCHO5MERAg4wuHOJKXhIKGA8P0/vEpHPE3AEQKO499E6Ey6dCbSlIaDNHf5AUHwA0ImduCpf+wMxAKCGal42sMRCgoI4F+lBByhI57mjnW7ilau8WxTmAfszXm8L7OsHxG5XkQ2isjGpqamYb1JV9LN/jrs4Yi/3Ex+e1tjlIQCACRdDxH/Uj/p+l+0klCArqRL2vOf6yGZq0Y3U0VRQE2FMaOu77lpKCKQ9jz2tcaKUyAmSUOzqt6pqqtVdfXMmUOO0u6lLBzo94XvqQc2k9+CmlK6U36ADwccv95X/fsA3Sm/HSno+M/1UD1aNQTD/3IaMxqG+2NEM1eu82tKi1MgxjcoNAALch7PzywbVX9x7mI89aOrp17mr7/cTH43nLeElKvEkmnqysO4qrieUlceJpZMk3KVvzh3MRXRIK6nuJ7n31QpDwcojwSJJdPUloWAzCW6HG1QdDK/zAYyKX5VmQkpGnTwFMrDhR1Fgn9lWxENcsN5S4pWrvHMfbQW+JiI3I3f0NymqgdG+016GpOt99HUdMHyem7Fb1vY1xrjhJlliAidiTT1FdFs76NT5lf36n10Ql3v3kfJtMf86kC291Fp0O99FE+5JF2lqmTg3keV1vtoxKz30fTqffQ/wAVAHXAI+CIQAlDVf890Sf1X4FL8LqkfUtUhuxUNt/eRMcaYwnsfFe1KQVXfPcTzCvx1sd7fGGPM8FmVqDHGmCwLCsYYY7Im3SQ7xhhjCpdMe9lxO4WwoGCMMVNIPOVmbh6JtIvrKaFA4ZVCFhSMMWYSS7sesZRLd9K/ecfYo9SCgjHGTDLxlEss6RJLpkmmRzcDlwUFY4yZBOIpl65Emq6EO+go+2NlQcEYYyaoZNqjM5GmK5EmNYzG4h6JlMvTu1tYt+Nwwa+xoGCMMROI5ymdyTQd8TSJ1PCzOaddjw27W3l0WyN/fKU5mzCyUBYUjDFmAuhOunQkUnQlXIabfkhV2Xawg4e3HuLx7U20daeyzwUc4YzjathT4LYsKBhjzDjxPKUjkaa9OzWi6qG2WIqHXzrEr188wO7mo3MsCHDqwmouPLGec5fWUVce4Z6/KmybFhSMMWaMpV2Ptu4UHfH0sLuQqipb9rez9oX9PLGjqddUsUvqyrh4xSwuXF7PzIrIiMpmQcEYY8ZIMu1xpDs5oiqiRMrlkZcauf/5Bl5p6souL4sEeMvyWVy+ajZL68sRObYZoywoGGNMkSXTHkdiSToTA8+ZMJCWriQPPN/A2hcO9GorOHF2BWteN5cLTpxJNDR6M0laUDDGmCJJux4tsSSdg0ygM5D9R7r53417+c3mg9kqoqAjXLi8nqtPm8vy2ZWjXVz/PYqyVWOMmcY8TznSnaKtOzXsaqI9zV389/rX+N32xuwczpXRIFe+bi5XnzqX2vKRtRUUyoKCMcaMovZ4itauJG4B04nmeq05xo/W7+F32xqz037WV0R45+r5XL5qzqhWEQ3GgoIxxoyC7qRLc1di2LmIDrbF+cFTu3nkpUPZK4N51SW896yFvOWk+mFlOB0NFhSMMeYYpFyP5s4kseTw2g3aYin+++k9rH1hf7bNYG51lD8/+zguOmkWAefYehHlioQsdbYxxhSVqnIkluLIMNsNkmmP+55r4L/X76Er6aegqCsP84FzFvHWk2cRHKUrAxGhPBKkqiREOGhBwRhjiqY76XK4MzGsUciqylOvNPO9J15h/5E4AOWRIO8+cwHXnDZv1NoMHBEqS0JUlYRGdLVhQcEYYwrkeUpzV5KOeGrolXPsa43xnUd3snFPKwCOwJpT5/H+c46jqiQ0KmULOg5VJSEqokGcY6h6sqBgjDEFGMnVQSLl8tM/vcbdG/Zm2w1ev6iGj15wPItqy0alXEHHoao0RGU0eMyjmcGCgjHGDEpVaY2lOBJLDut1z+5p5R9/u4MDbX5V0azKCB978wm88YS6USlX0HGoLgtRERmdYJDd7qhtyRhjppiU69HYkRjWvAYd8RT//sQufr35IOCPQn7X6xfw3rMWjkq7QcARqkvCVJaMbjDoYUHBGGPy6EqkaepIDCuL6VOvHOaffvsyLV3+VcXJcyv51CXLRqWqSESozjQgH0ubwVAsKBhjTB8tXclhVRd1JtJ893c7eWjLIQBKQgH+4k2LWXPqXJxR+DVfHg0yozQ8at1VB2NBwRhjMlxPaeyI050svLroudda+fpvttPYkQDgtIXVfPqtJzK7MnrM5YmGAtSWh4kExybFBVhQMMYYwB9Udqg9XnDvorTr8V9P7ebuP+1FgUjQ4frzlozK1UEo4FBTFqY8MvanaAsKxphprzvpcqg9XnD7QUNrN1958CW2H+wA/LkNPnfZchbMKD2mcvS0G1SXhorSiFyIogYFEbkU+DYQAL6vql/r8/xC4IdAdWadm1T1wWKWyRhjcrXHUzR3JgtOVfG7bY188+EddKdcBLjuzAV86A2Ljrm+vzQcpLY8POYJ8PoqWlAQkQDwXeBiYB+wQUTWqurWnNU+D/yfqn5PRFYADwKLilUmY4zJNZwG5WTa43tPvMIDz+8HoLYszGcvW87px9UcUxmCjkNteZiycagqyqeYpTgT2KmquwBE5G5gDZAbFBTomT6oCthfxPIYYwzgD0hr6kwUPCPawfY4f/+LrdnqojOOq+Fzly+npjR8TOWoLAkxozRc1C6mw1XMoDAP2JvzeB9wVp91vgQ8LCIfB8qAi4pYHmOMwfOUQ8PoYfTsa63c+outtMfTCPCBNxzHe8867phSW4eDDnXlkTGbOGc4xvt65d3AD1T1H0XkHODHIrJSVXs1/4vI9cD1AAsXLhyHYhpjpgLXUw60dRc0EY6q8rNn9nHHul146k+J+YUrVnDGMVQXiQg1pf4AtPFqSB5KMYNCA7Ag5/H8zLJcHwEuBVDVP4pIFKgDGnNXUtU7gTsBVq9ePbw57owxBj9lxcG2wrqcJtMe33x4O4+85J+KTqgv59Y1Jx/T2INIKEB9RWTcG5KHUsygsAFYKiKL8YPBdcB7+qzzGvAW4AcichIQBZqKWCZjzDSUSLscakuQ9oYOCC1dSb7wwGZeOuC3H1x0Uj2fungZkRFW9fR0M60pO7b2h7FStKCgqmkR+RjwEH5307tUdYuI3ApsVNW1wKeA/xCRv8VvdP6gDmcKI2OMGUI85Y9BcL2hTy2vNHZy8/2baexIIMAN5y/hHWfMH3FVTyjgMLNiYrYdDKSobQqZMQcP9ll2S879rcAbi1kGY8z0FU+5HGwrbFDa+l3N3PrLrcRTHiWhADe/bTlvOH7kaa4roiHqysMTtu1gIOPd0GyMMUUxnFHKv9y0n2898jKeQn1FhNuuWcnxM8tH9L4BR5hZEaE0PDlPr5Oz1MYYM4jupMvB9viQo5RVlbue3M1Pnn4NgKX15Xz12lXMGGH9f2k4yMyKyDF1Vx1vFhSMMVNKTxvCUAEh7Xp88+EdPLzVT3d91uIZ3HLFCkrCw6//FxFqy8NURkdnvuXxZEHBGDNlFNqGEE+53PrLrazf1QLA21bN4W8uWjqiX/jhoEN9RZRwcGJ3NS2UBQVjzJSQSBcWEDriKW6+bzOb97cD8P5zjuMD5xw3ogbh8miQmeWRSdeYPBgLCsaYSS/tehxqG3rqzJauJJ/52SZ2He5CgI9deALXnDZv2O83laqL+rKgYIyZ1PzUFfEhB6Y1tse58Web2NfaTdARbrpsORcurx/2+wUdh/rKyTX2YDgsKBhjJi1VLWi2tIYj3dx4zwscak8QDjp86coVnL2kdtjvFwkFmFURGZO5kseLBQVjzKTV2JEgnho82+me5i5uvGcTzV1JoiGH265eyWkLh5/Ubiq2H+RjQcEYMykd7kzQlRh8PoTdzV186v9eoDWWoiwS4GvXruLkuVXDeh8RYUZZmKqSqdd+kI8FBWPMpHMklqS9OzXoOq8e7uLGe/yAUBENcvvbT2HZrIphvU/AEeoroiMauzBZWVAwxkwqHfEULV2DT6H56mH/CuFId4rKTEBYOsyAEAo4zK6KTvhU16PNgoIxZtLoTroc7hw8IPRUGfUEhH98x+s4vn54eYyioQCzKqOTOl3FSFlQMMZMCsm0N2T6ir0tMW68Z9PRgPDO1w07sV15xM9fNNUblAdiQcEYM+G5ng6Z8bThSDefuucFWrqSlEf8KqPhBoTq0vCIk+FNFRYUjDETWiFjEQ61x/nU/73A4c4kZeEA33j7qmG3IdSWRagqnR49jAZjQcEYM6E1dQ4+FqGlK8mnf7aJxo4EJaEAX/uzVSyfXVnw9kX8+Q/KI3Y6BAsKxpgJrC2WojM+8FiEjniKz/zcT10RDjrcds3KYY1DcESYVTm9upwOxYKCMWZC6k66NHclBn3+s/e+yK6mLgKO8KUrV3DqguqCtx9w/IAwVXMYjZQFBWPMhJNyPRo74gM+n0x73PLAZrYe6ECAz122fFi5jIKOPwZhqsyBMJosKBhjJhQv09PI9fL3NHI95au/3sYzrx0B4JMXL+PNw8h2Ol0HpRXKgoIxZkI53Jkgmc7f00hV+c5jL/PEjiYA/vJNi3nbKXMK3nY46DC7Mjqls5weKwsKxpgJ40gsSecgSe5++NQefvHCAQDeccZ8rnv9goK3HQkFmD1NRykPhwUFY8yE0J10B81ptPaF/fxo/R4ALlkxixvOX1LwqOOScIBZFVEcCwhDsqBgjBl3QzUs/+Hlw3zn0ZcBOHPxDG68ZBlOgQFhuqetGC4LCsaYcTVUw/Lmhja+8uBLeAonzq7gi1euKLhNoLIkRF15ZDSLO+VZUDDGjKumQRqWX2uJ8fn7N5NMe8yrLuEfrllJSYHjCmpKw9RM8zxGI2FBwRgzblq7kgPOntbSleSz975IezxNTWmIr/3ZKmpKCzvJ15ZHps1MaaPNgoIxZlx0JdK0xvI3LHenXG6+fzMH2uJEM+kr5lWXDLlNEaGuPExF1ALCSFlQMMaMuWTao6kjfwoL11Nu+9VLbD/YgSPw+StOKijBnYgwqzJCadhOa8fC9p4xZkx5Q8yN8L3HX+GpV5oB+NibT+ANx9cNuU1LbDd6ijqsT0QuFZHtIrJTRG4aYJ13ishWEdkiIj8tZnmMMeOvsSMx4NwI9z7bwL3PNQD+4LSrT5s35PYcEWZXWUAYLUW7UhCRAPBd4GJgH7BBRNaq6tacdZYCnwXeqKqtIlJ4AhNjzKTT0pUklszfsLx+VzP/9vhOAN60tI4bzl8y5PYs0+noK+aVwpnATlXdpapJ4G5gTZ91/hL4rqq2AqhqYxHLY4wZR12JNEcGaFje2djJrb/cmh2L8NnLlg85OK3nCsECwugqKCiISKmIfEFE/iPzeKmIXDHEy+YBe3Me78ssy7UMWCYiT4rIehG5dID3v15ENorIxqampkKKbIyZQAZrWD7cmeBz971IPOVRXxHhtqtXDnmi7wkIkaAFhNFW6JXCfwEJ4JzM4wbgK6Pw/kFgKXAB8G7gP0Sk3ywZqnqnqq5W1dUzZ84chbc1xowVVaWxI3/Dcjzl8vn7N3O4M0lpOMBXr13FjCEGnIlYlVExFRoUjlfVbwApAFWNAUMlEmkAclMYzs8sy7UPWKuqKVV9FdiBHySMMVPEQCOWPVW+9utt7DjUiSNwyxUrWFxXNui2RIT6iog1KhdRoUEhKSIlgAKIyPH4Vw6D2QAsFZHFIhIGrgPW9lnnfvyrBESkDr86aVeBZTLGTHDt8YHnWP6vJ3ez7uXDAPx/F5zAmYtnDLqtnoBQFrGe9MVU6N79IvAbYIGI/AR4I/DBwV6gqmkR+RjwEBAA7lLVLSJyK7BRVddmnrtERLYCLvBpVW0e2UcxxkwkibRLc2f+huXfbj3ET55+DYCrXjeXa06bO+i2bGDa2BEdYABJvxVFaoGz8auN1qvq4WIWbCCrV6/WjRs3jsdbG2MK5HlKw5HuvOMRtuxv45P/9wIpVzljYTVfvXbVoFlPRYTZNjDtmInIM6q6eqj1htMldR7+L/4wcJ6IXDvSwhljprbDnfkHqB1qj3PLA1tIucr8mhJuGSINtgWEsVfQtZiI3AWcAmwBev7TCtxbpHIZYyap9ngq75Sa3SmXL9y/hdZYivJIkNuuXjlo4rqeKiMLCGOr0Aq6s1V1RVFLYoyZ9AZqR+jpabSzye9p9MUrV7BgRumA27E2hPFTaPXRH0XEgoIxZkCepzS2J8jXTvnDp3bz+0xPo49feAJnHFcz4HZ6ehlZQBgfhe71H+EHhoP4XVEFUFU9pWglM8ZMKk0DtCM8vr2RH68/2tNozamDJ7mbad1Ox1Whe/4/gT8HXuRom4IxxgDQFkvlnUFtx6EOvv6b7QCcuqCKj735+EG3U18ZpdwCwrgqdO83ZcYVGGNML/GUS0ueRHctXUm+cP8WEmmPOVVRvnjFyYP2NJpZEbGAMAEU+h94LjPXwS/IGcmsqtb7yJhpLO16edsRkmmPWx7YQlNngpJQgK9cvZKq0oF7GtWWR2wKzQmi0KBQgh8MLslZZl1SjZnG/ER3CdKe12/5tx55ma0H2hHgc5cvHzSn0YyyMFUlFhAmioKCgqp+qNgFMcZMLs1dSeIpt9/ye59r4DdbDgLw4XMX8cYTBp5Os7o0THXp4FlRzdgqdD6F+SJyn4g0Zm4/F5H5xS6cMWZi6oinaO9O9Vu+cXcL33v8FQDefOJM3nPmwgG3UVkSGjJNthl7w5lPYS0wN3P7RWaZMWaaSaRdDucZoLavNcaXf/USnsLS+nI+/dYTkQFmTyuPBKkrjxS7qGYECg0KM1X1v1Q1nbn9ALDZboyZZgYaoNaVSPOF+7fQEU9TUxriy2tOHnASnNJwkJkVFhAmqkKDQrOIvE9EApnb+wBLcW3MNNPY0X+Amusptz34EntaYgQd4e+vOpn6ymje10dDAWZVRga8gjDjr9Cg8GHgncBB4ADwdsAan42ZRo7EksSS/Qeo3fXkq6zf1QLA31y0lJXzqvK+PhIKMLsyagFhgiu099Ee4Koil8UYM0F1J11auvq3Izz6UiP/86e9AFx72jwuXzUn7+sjoQBzKqM4jgWEia7Q3kc/FJHqnMc1mXTaxpgpLu16NHbE+y3fcaiD2x/2U1icvrCaj16QP4WFBYTJpdDqo1NU9UjPA1VtBU4rTpGMMRNJY0cC1+vdsNzSleTz928mmfaYWx3lC1esIJDnpB8OOsy2gDCpFBoUHBHJ5roVkRkUPhraGDNJteQZoNaTwuJwZ5KSUIAvr1mZd0RyKOAHhHzBwkxchZ7Y/xE/dfY9mcfvAP6hOEUyxkwEsWSaI30S3akq//zIjmwKi5vflj+FRcARZlVGB02AZyamQhuafyQiG4ELM4uuVdWtxSuWMWY8pVyPpo5Ev+U/e7aBh7YcAvwUFm84vn8KC3/WtCjhoAWEyajQOZp/rKp/DmzNs8wYM4WoKofa4/3aETbsbuGOJ/wUFhcurx8whUV9RWTAgWtm4is0lJ+c+0BEAsAZo18cY8x4a+pMkEz3HqD2WkuMW3+5FU9h2axyPn3JsrzjDeps1rRJb9CgICKfFZEO4BQRaReRjszjRuCBMSmhMWbMtMdTdMZ7D1DriKf4/P2b6Uq4zCgL8+U1K4nkuRKoKQ1TaXMiTHqDBgVV/aqqVgC3q2qlqlZkbrWq+tkxKqMxZgzEUy7NfRLduZ5y6y9fYl9rN6GAcOtVJ+fNW1QRDVFjGU+nhEKv834tIuf1Xaiq60a5PMaYceB6SlNH/0R3//b4KzyzpxWAGy85kRVzK/u91hLcTS2FBoVP59yPAmcCz3C0N5IxZhJrypPo7hcv7Oe+5xoAePeZC7h4xax+rwsHHeotIEwphXZJvTL3sYgsAL5VlBIZY8ZUa1f/RHfPvdbKdx7bCcA5S2r5yLmL+70u6Nho5alopB2J9wEnjWZBjDFjL5ZM09pngNq+1hhf+sVWXE9ZXFfGzW9bjtOnp5EjwqyqiA1Om4IKHafwL0BPZaODn/fo2WIVyhhTfMm0R2N77wFqHfEUn7tvMx3xNNUlIW67eiWl4d6nCRGhvjJCJGhjEaaiQsP8VmBH5rYe+Iyqvm+oF4nIpSKyXUR2ishNg6z3ZyKiIrK6wPIYY46B5/kD1LychuW063HrL7Ye7Wm05mRmV/WfLKeuPNwvUJipY9D/rIgE8XMcfRh4LbN4IXCXiPxJVfvP3H30tQHgu8DF+NVNG0Rkbd/0GCJSAXwCeHrEn8IYMyx9Z1BTVf71d6/wzGt+MuQbLzkx72Q5tWURKmwswpQ21JXC7cAMYLGqnq6qpwNLgGrgm0O89kxgp6ruUtUkcDewJs96Xwa+DvRP2G6MGXXNnYl+Dcs/f7aBtS/sB+C9Zy3M29OoqiREVakFhKluqKBwBfCXqtrRs0BV24GPApcP8dp5wN6cx/syy7JE5HRggar+quASG2NG7EgsSVt37wv8p145zPce93Manb9sJh9646J+ryuPBqktt66n08FQQUG172gWf6HL0YbnERERB/gn4FMFrHu9iGwUkY1NTU3H8rbGTFvt8VS/KTVfPtTBV371EgqcNKeCmy49sV9Po7JIkPqK/m0LZmoaKihsFZH3910oIu8Dtg3x2gZgQc7j+ZllPSqAlcDjIrIbOBtYm6+xWVXvVNXVqrp65syZQ7ytMaavrkSaw31SYTd1JLj5/s3EUx6zKiN5cxqVhAM2OG2aGaoLwV8D94rIh/FHMAOsBkqAa4Z47QZgqYgsxg8G1wHv6XlSVduAbDJ2EXkcuFFVNw7nAxhjBhdPuTT2CQhdiTSfve9FDncmKQ0H+IdrVjGjT+6iSCjArIpo3myoZuoaNCioagNwlohcyNH02Q+q6qNDbVhV0yLyMeAhIADcpapbRORWYKOqrj3GshtjhpBMexxqj/fKaZR2PW795VZ2NXURcIS/v+rkfrOn2dzK01ehaS4eAx4b7sZV9UHgwT7Lbhlg3QuGu31jzMBcr/9kOarKtx/dyYbdfpK7T160lDOOq+n1ulDAYU5Vic2tPE3ZGHVjpiDPUw62x/slufvpn17jVy8eAOB9Zy/kslVzej0fCjjMrbaAMJ1ZUDBmilFVDnXESaTcXssf3nKQ//zDbgAuOqmeD71hUa/ng47DnKqoBYRpzoKCMVOIP79ygu5k74CwYXcLtz+8A4DTF1bz6bee2KsB2RLcmR52BBgzhTR19B+t/PKhDr601s96umRmGV+66mRCOSd/EWFWZdQS3BnAgoIxU0ZjR5zORO+A0HCkm5vufZHulEt9RYSvXbuK8kjv/iUzKyKUhC0gGJ8FBWOmgKaOBJ3x3gGhpSvJ3/18E62xFBXRIF+9dhV1fVJV1JZH+gUJM71ZUDBmkmvqSNAR753PqCuR5qZ7X2T/kTiRoMNtV6/sNxahpjRMVYkluDO9WVAwZhLLFxCSaY9b1m5hZ2MnjsAtV6zolwa7qiRETZ8RzMaABQVjJq18AcH1lC//aivP5cyLcM7xtb3WsYynZjBWmWjMJNTYEe/XhuCpcvtD23lyZzMAf3X+Ei5dObvXOpbx1AzFgoIxk0xje/9eRqrKd3/3Cg9vPQT4E+W8c/WCXutYxlNTCAsKxkwSqkpjR4KuPgEB4K4nd3Pfc35m+jWnzuXDfSbKsYynplAWFIyZBHpGKvcdmAbw3+v38JOn/SnULzqpno9feEKvk79lPDXDYQ3NxkxwgwWEezbu5a4ndwNw3tI6/u7S5b1mTgsHLeOpGR67UjBmAvM8P7ld31xGAA8838D3ntgFwNlLZnDz207qdfK3gGBGwoKCMRNU2vU42B4nmfb6Pbf2hf18+9GdAJxxXA1furJ3PiObE8GMlAUFYyagZNrjYFuctJc/IHzrkZcBOHVBNV9eczLh4NGAYCmwzbGwoGDMBBNPuf1mTOvROyBUcds1K4mGjiazCzjC7KqopcA2I2ZBwZgJpDvpBwRP+weEA8RIAAAaDklEQVSE+55r4F8e86uM/ICwipKcgOBkUmDnXjUYM1wWFIyZILoSaRo7EmiegHD3hr3cuc5vVM4XEESE+spIr6sGY0bCgoIxE0B7PMXhjkS/5arKj9fv4QdP7QFg9XE13Lrm5H4n/5kVEUrD9nU2x86OImPGWUtXkiOxZL/lqsp//P5V7t6wF4A3HF/LLVes6Fc9ZHMimNFkR5Ix40RV/clx8qStcD3lW4+8zK9ePADA+ctmcvPly/s1IFfbnAhmlFlQMGYcpF2PQx0JEqn+g9KSaY+v/nobT+xoAuDylbP524uX9etiWlUSYobNiWBGmQUFY8ZYLJmmqSORt8tpLJnmi2u38syeVgDetXo+15+3pF8iu6qSkM2JYIrCgoIxY2ig9oOe526690V2NnYC8JdvWsy7z1zYb71KCwimiCwoGDMGXE9pHCCHEcBrLTFu+vmLHGyP4wh88uJlXL5qTr/1KktC1FlAMEVkQcGYIounXBrbE3lTVgC8uK+NLzywmfZ4mmjQ4ZYrV3D2ktp+61WXhq0NwRSdBQVjiqitO0VLVzLvgDSA3249xDcf3k7KVapLQvzDtStZPruy33o1pWFqLCCYMWBBwZgi8DylqTP/LGngd0f9wVO7+fF6f3KchTNK+YdrVjK3uqTfurVlEapKrdupGRsWFIwZZfGUS1NHgpSbv7qoO+Vy+2+283imy+kZC6v54pUnUx7t/3W0gGDGWlEzZ4nIpSKyXUR2ishNeZ7/pIhsFZFNIvKoiBxXzPIYU2xHYkkOtMUHDAgH2+L8v/95LhsQrjhlDl+9dlX+gFBuAcGMvaJdKYhIAPgucDGwD9ggImtVdWvOas8Bq1U1JiIfBb4BvKtYZTKmWIbqXQTw3Gut/P0vttIeT+MI/PWbT+DqU+f2G4MAUFcRoTJqAcGMvWJWH50J7FTVXQAicjewBsgGBVX9Xc7664H3FbE8xhTFYIPRwG8/+N+N+/j+73fhKVRGg3zxyhWctrCm37oiQn1FhDLLZWTGSTGPvHnA3pzH+4CzBln/I8Cv8z0hItcD1wMsXNh/MI8x48H1lObO/LmLenQm0nz9N9t4cmczACfMLOfWNSczuyrab92A48+HYOmvzXiaED9HROR9wGrg/HzPq+qdwJ0Aq1evzv9zzJgx1BH3u5oOdHUAsLOxk7//xVYajnQDcOnJs/nEW04gkuekHwo4NkGOmRCKGRQagAU5j+dnlvUiIhcBNwPnq2r/hPLGTCCJtEtzZ5J4nkR2PVSV+55r4I51u0i5SiggfOItS/OOUAYoDQepr4jg2JzKZgIoZlDYACwVkcX4weA64D25K4jIacAdwKWq2ljEshhzTDxPaYklae9ODbpeWyzFNx7azh93+dVF86pL+MIVJ7FsVkXe9W1QmploihYUVDUtIh8DHgICwF2qukVEbgU2qupa4HagHLgn0wPjNVW9qlhlMmYk2rpTHIkNXlUEsGF3C994aDvNnX7Cu4tXzOITbzkh74xoAUdstjQzIRX1iFTVB4EH+yy7Jef+RcV8f2OORXfSpbkrQTKdf8xBj3jK5Y51u3jg+f0AlIQCfOKipVyyYlbe9SOhALMqIv0mzDFmIrCfKcb0kUi7tHaliCUH7lXUY3NDG994aDv7Wv3G5JVzK7npsuV501XA0Ylx8o1NMGYisKBgTEba9WiJJemMDx0MulMu//n7V7nvuQYUCDrCB9+wiHe9fkG/GdIAHPGri2z8gZno7Ag1057rKa2xJB3x9IDZTHM9s6eVf/rtDg60xQE4ob6cz7z1RE6oL8+7fiQUoL4iQsiqi8wkYEHBTFuup7R1p2jvTuEVEAxaupL82+Ov8Ng2v6NcKCB84JxFvHP1/AHbB6y6yEw2FhTMtDPcYOB6yi837ef7f3iVroQ/PmHVvEr+9uJlLKoty/uaoOMwsyJCSdhGJ5vJxYKCmTa8TDBoKzAYADy/9wj/+rud7GrqAvy8RTect4S3rpyNM8Cv//JIkNrySN62BWMmOgsKZspTVdq70xzpHnqsQY8Dbd3csW4X63YcBkCAS1fO5vo3LRkwnXXQcagtD1tjspnU7Og1U9ZIgkFbd4qfPL2H+5/bTzrzmhVzKvn4hSdw4uz8o5IBKqIhasvClqrCTHoWFMyU43lKe9yvJio0GHSnXO57toH/2fBatt2gviLCR85dzEUn1Q/YUBwOOtSVRyyzqZkyLCiYKWO4DcgAiZTL2k0HuPtPr9Ea8/MalUUCvPfMhVxz2ry8GU3Bn/egpjREVUnIehaZKcWCgpn0EmmXtu4UXQm3oHEG4F8Z/HLTAf5v495srqJw0GHN6+bynrMWUlUy8Kxn5ZEgM8rClqbCTEkWFMykpKp0JtK0x9MkBklj3VdnPM0DLzTws2caaMtkPA06whWnzOE9Zy2krjwy4GsjoQC1ZWGrKjJTmgUFM6kk0x4d8RSdiXTB7QUAB9vi/OzZffz6xYN0Z4JIKCBcvmoO171+AbMq+8+E1iMUcKgpC1NuvYrMNGBHuZnwPE/pTKbpGOZVgaqyaV8b9z3fwB9ePkxPDCkJBbjqdXN4x+oFzBhkLoOg41BdFqIiErR2AzNtWFAwE5KqEku6dCXSdCULbysAiCXTPPpSI/c/v59XD3dll9eVh7n29PlcsWoO5dGBD/1QwKGq1IKBmZ4sKJgJQ1WJpzw6E2liyeFVD6kq2w528KtNB3hseyPx1NE5EFbMqeDq0+Zx/rKZgyalCwUcqktDlFswMNOYBQUzrlSV7pRLV8IddiAAaOpI8MhLh3h46yH2NMeyy8NBhzefOJOrT5036KAzgJJwgKqSkM2CZgwWFMw4cD0/EMQSaWJJt+AxBT3aulP8/uXDPL69kedeO0Luq4+fWcbbVs3hLSfVUxEduFupiFAeCVJZEiQStN5ExvSwoGDGRCLtEk96dCXTJNLesNoIANpiKZ565TBPvHyYZ/a09rqiqCoJ8Zbl9Vxy8iyW1pcPWvUTCjhURkOUR4OWsM6YPCwomKLouRroTvq3tDf4PMf5NLR2s/7VZp7ceZhN+9rIrVmKhhzecHwdFy6fyZmLZgw6kExEKA0HqIyGLJW1MUOwoGBGhecp8XQmCKTcISe7zyeRcnmxoY0Nu1tZv6uZvZl5j3tEQw5nLa7l/GUzOXvJjCEHkYWDDhURuyowZjgsKJgRcT0lnnL9W9ojOYIqIddTdjZ28tzeIzy7p5VNDW39gklNaYizFtdy7tJaVh83g3Bw8NQS4aBDWThIaSRgbQXGjIAFBTMkVSXpesRTHom0SyLlkXKHfyWQcj12HOrgxX1tvNjQzqaGI9mMpD0EWDa7gjMX1XD2klpOnF0x4GQ2PaKhAGXhICXhwJBBwxgzOAsKppeeAJDI/PpPjPAqQFVp7Eiw/WAHWw+0s3V/OzsaO/NWK82pinLagmpWL6rhtIU1gyajA3+kcTTsUBoOUhIKWNWQMaPIgsI05npKMu2RdL1ef0cSAA60xdnZ1MkrjZ3sONTJjkMd2VTUfc2vKWHl3CpOmV/FqQurmT1I3iEAR4SScIBoMEA07Fi1kDFFZEFhGsg9+acyt2TaG/ZAMVXlSHeKPc0xdh/u4tXmLl5t8v/2rQbqEQ46nDirnBVzKjlpbiUr51YNmm8IIOAI0ZAFAWPGgwWFKaKn2iftqn/Sdz1SrpJ2h3/y7065HDjSTcOROPtaY+xr7WZfa4w9zTHa4+kBXxcOOiypK2NpfTknzq7gxFkVHFdbOmR30XDQIRxwiIYcoqHAoKkojDHFZUFhEkm7HmnPP+m7npLKBIC0q8MaB5BMezR2xGlsT3CoI8HBtm4OtMU51B5n/5E4zV3JIbcxqzLC4royltSVsbiunBPqy5hfUzpo/X7AORoAQkGHSOa+5RkyZuKwoDABqCqup6S9Pn8zQaBn2VB1/WnX40h3ipauJC1dSZo7kzR3JWjuTNLUmeBwR5LDnQmOdOev6+8rEnSYX1PC/JpS5teUsHBGKcfVlrJgRiklg0xTGQoIoYBDKOAQDIgfBAKONQgbMwlYUCgSz1PczMne6/nrQdrzsstzb325nj+zWEc8RUc8TXs8RXt3Ojshfc/tSCxFWyxFayw5aNVOPo5AbVmE2VURZleVMLvS/zuvOsrc6hJqy8L9fsUHHCHgCEHHP+EHHSEYcAg6Yid+Y6aAogYFEbkU+DYQAL6vql/r83wE+BFwBtAMvEtVdxezTMOhqnhK9sSuCq5m7ntkT+6qfgBIuR6xhEtX0qU7mSae8kf3xjKjfLuTmftJl65kOrNums5Emq6e+/F0dg6BYxEKCLVlEWrLw9SWhakrj1BXEWFmeZi6igizKqPUZeYZFvFP7o4jBESyJ37/5N/7vlX1GDO1FS0oiEgA+C5wMbAP2CAia1V1a85qHwFaVfUEEbkO+DrwrkK233PCTnuZxtWeLpWuRyqtJF2XlOv//eeHtvOHV1pQ/MFRpy+s5l2vX0DSVZJpl0ROf/yeXjq5jxNpN9tdM5H2SKT8+/G0n84hd1RvsQQdobIkRGU0SFVJKHurLg1RXRqmpjTMjLJQ9uRfGQ0ScBwcERzxf+FL5oT/5MtNfPuRl9nXGmNBTSl/df7xXLC8fthlenxbI3es28WW/W10JV08T6mIBvmLcxfz/y5axuPbGvn6b7axKzPRzeLaUm667CQA7li3i72tMSoiQfa2dNGZ9PedI7DmdXM42Bbnj6+2jt4ONBNKyIFZlVEOtMVxh+gH4QDRcIBYzg8lEZhREqS+soTGjjgpVwkHHWaWR1BV9rTE6M6ZU6OuLMQ333HqgMd5z7G8tzVGeTiAiNCRSLOgppQbzlvCBcvr+c4jO/j3dbuIJV1EYF5lhK9ccwpw9HjOXX+g7VdEgv4c40k3e7+pM5H9DEvrKzhnyQz+uKtl0G0Wiwy3T3rBGxY5B/iSqr418/izAKr61Zx1Hsqs80cRCQIHgZk6SKGic5bq/A99m/Qwe9SMp4BASThIaTiQuflpGCoiQcoiQcozt4rMSb88EqS6JERlSYiq0hDVJSHKIn7+HkEQIXuyl5y/hXh8WyO3rN1CKCCUhAJ0p/zgeetVJw/roOvZTmc8RUvOeISeslx1ymz+sLOZ1liKnholT6Ek08OosiRE2vXY09I9wDsYU5iA438fvMw5YaAgUxpy+Lf3npH3hN3znUi7Hg1H4gDMq44SDDikXOWMhVU88MIB+p52okGhIup/Vwf6Pg20/ZrSIK2xNJ6nSOYzoFARDdIWT/tX9eWREX9H+xKRZ1R19VDrFbP6aB6wN+fxPuCsgdZR1bSItAG1wOGBNqowagGhtix8tDdMQAgH/TQJkaDjLw86RIMBIkGHSMjvLx8J+ie1knCAklCAaMihJBSgNBKkNOSf9MsifsqF0kwgiAQnTg+bO9btIhSQ7IQypeEgsWSaO9btGtYB17OdI90pBECgJ5Q7Ams3HSSUqXbqSVMhntKV9A/w2VUl7GrqHOVPZ6Yjz4NwyCHuuTDIqSGW8vIe57nfiV1NnQREQOBwZ5IlM8uJJdOs3XQQT8ke6+Af7/G0QjzN7KoSIP/3qd/2Hf/k39yVIhRw/CrozGfwMmOBAo7QEU8zsyI64u/oSE2KhmYRuR64HmDW/EX86MNnZuu6exo5A5mGzmBACGUaQUMBh9ff9kifbR09eT3zhYvH+qOMu72tMar7pJEoCQXY1xob4BWDbyf7RclQ9YNCylUE7TVGQYRslR9AcgT5k4zpS3PuDPVzMd9xnvudSLpeNij0HJ8locDRziB5ftv17Q7e9/vUb/uZS2dPM+ejnNf2fEdC0vv7MZLv6EgVMyg0AAtyHs/PLMu3zr5M9VEVfoNzL6p6J3AnwOrVq/W8ZTMLLkTAEVxP6ftDfbr2kllQU0pjR7zX1JPdKZf5NaUj2o6TE2Th6EHd0zCtmQMfjgaMoOMHinDAIeUeW4O6MZJzR4YIDPmO89zvRDjgkM7UP4UzP2i6U272PJJtmMzRczz36Pt96rf9zHZ6vjv5flR5evT9822zmIo5dHQDsFREFotIGLgOWNtnnbXABzL33w48Nlh7wkhcdcpswN/ZPbfc5dPNDectIeUqsWQaVf9vylVuOG/JiLZTXRJC6R0YPPX3b0U0mOly6/k3VcrCASqi/uVwXfng6S6MKYTj+L/WHWCw33qlISfvcZ77nagrD2d7FdaVh7Pfj6tOmU2m1qfXecRvUwgO+n3qt/1Md/XaspD/gzXnM3ieZq/AK6LBY/qOjlTRGpoBRORy4Fv4XVLvUtXbRORWYKOqrhWRKPBj4DSgBbhOVXcNts3Vq1frxo0bh1WOv737WdZuOojrKQHHbwT95+tOH9Fnmgp6ekLsa40x/xh6Nhxr76N9rTHKrffRtDSavY+aOuIkR6n30b7WGGWZ3kediXSv78dQvY8G+z7lbr880+OoK+lm7x/uTGQ/Q27vo2P9juYqtKG5qEGhGEYSFIwxZrorNChY5jFjjDFZFhSMMcZkWVAwxhiTZUHBGGNMlgUFY4wxWZOu95GINAF7RvjyOgZJoTHOJmrZrFzDY+UavolatqlWruNUdciRv5MuKBwLEdlYSJes8TBRy2blGh4r1/BN1LJN13JZ9ZExxpgsCwrGGGOypltQuHO8CzCIiVo2K9fwWLmGb6KWbVqWa1q1KRhjjBncdLtSMMYYM4hpExRE5FIR2S4iO0XkpnEsxwIR+Z2IbBWRLSLyiczyL4lIg4g8n7ldPg5l2y0iL2bef2Nm2QwR+a2IvJz5WzPGZToxZ588LyLtIvI347W/ROQuEWkUkc05y/LuI/F9J3PMbRKRoqXmHaBct4vItsx73yci1Znli0SkO2ff/fsYl2vA/52IfDazv7aLyFvHuFz/m1Om3SLyfGb5WO6vgc4PY3eMqeqUv+Gn7n4FWAKEgReAFeNUljnA6Zn7FcAOYAXwJeDGcd5Pu4G6Psu+AdyUuX8T8PVx/j8eBI4br/0FnAecDmweah8BlwO/xp9H5Wzg6TEu1yVAMHP/6znlWpS73jjsr7z/u8z34AUgAizOfGcDY1WuPs//I3DLOOyvgc4PY3aMTZcrhTOBnaq6S1WTwN3AmvEoiKoeUNVnM/c7gJfw56qeqNYAP8zc/yFw9TiW5S3AK6o60sGLx0xV1+HP/ZFroH20BviR+tYD1SIyZ6zKpaoPq2o683A9/uyHY2qA/TWQNcDdqppQ1VeBnfjf3TEtl4gI8E7gf4rx3oMZ5PwwZsfYdAkK84C9OY/3MQFOxCKyCH+Coacziz6WuQS8a6yraTIUeFhEnhF/XmyAWap6IHP/IDBrHMrV4zp6f1HHe3/1GGgfTaTj7sP4vyh7LBaR50TkCRF50ziUJ9//bqLsrzcBh1T15ZxlY76/+pwfxuwYmy5BYcIRkXLg58DfqGo78D3geOBU4AD+5etYO1dVTwcuA/5aRM7LfVL969Vx6a4m/pSuVwH3ZBZNhP3Vz3juo4GIyM1AGvhJZtEBYKGqngZ8EvipiFSOYZEm5P8ux7vp/eNjzPdXnvNDVrGPsekSFBqABTmP52eWjQsRCeH/w3+iqvcCqOohVXVV1QP+gyJdNg9GVRsyfxuB+zJlONRzOZr52zjW5cq4DHhWVQ9lyjju+yvHQPto3I87EfkgcAXw3szJhEz1THPm/jP4dffLxqpMg/zvJsL+CgLXAv/bs2ys91e+8wNjeIxNl6CwAVgqIoszvzivA9aOR0Ey9ZX/Cbykqv+Uszy3HvAaYHPf1xa5XGUiUtFzH7+RcjP+fvpAZrUPAA+MZbly9Pr1Nt77q4+B9tFa4P2ZHiJnA205VQBFJyKXAp8BrlLVWM7ymSISyNxfAiwFBp0bfZTLNdD/bi1wnYhERGRxplx/GqtyZVwEbFPVfT0LxnJ/DXR+YCyPsbFoUZ8IN/xW+h34Uf7mcSzHufiXfpuA5zO3y4EfAy9mlq8F5oxxuZbg9/x4AdjSs4+AWuBR4GXgEWDGOOyzMqAZqMpZNi77Cz8wHQBS+PW3HxloH+H3CPlu5ph7EVg9xuXaiV/f3HOc/Xtm3T/L/I+fB54Frhzjcg34vwNuzuyv7cBlY1muzPIfAH/VZ92x3F8DnR/G7BizEc3GGGOypkv1kTHGmAJYUDDGGJNlQcEYY0yWBQVjjDFZFhSMMcZkWVAw056IzBKRn4rIrkyKjz+KyDUicoGI/HK8y2fMWLKgYKa1zGCh+4F1qrpEVc/AH9w45snjjJkILCiY6e5CIKmq2Rz5qrpHVf8ld6XMHAA35jzenElYhoi8P5Pc7QUR+XFm2SIReSyz/FERWZhZ/o7Ma18QkXWZZQHx5z7YkFn/hqJ/amMGEBzvAhgzzk7GH6U6IiJyMvB54A2qelhEZmSe+hfgh6r6QxH5MPAd/HTHtwBvVdUGyUx6gz/Kt01VXy8iEeBJEXlY/fTRxowpu1IwJoeIfDfzK35DgS+5ELhHVQ8DqGpPjv5zgJ9m7v8YP30BwJPAD0TkL/EnDQI/z9T7xZ/p62n8lAZLj+2TGDMydqVgprst+LltAFDVvxaROmBjn/XS9P4RFR3Jm6nqX4nIWcDbgGdE5Az8/DUfV9WHRrJNY0aTXSmY6e4xICoiH81ZVppnvd340zeSmQd3cc7r3yEitZnneqqPnsJvsAZ4L/D7zPPHq+rTqnoL0ISf9vgh4KOZlMmIyLJMplpjxpxdKZhpTVVVRK4G/llEPoN/ou4C/q7Pqj/Hr+LZgl/FsyPz+i0ichvwhIi4wHPAB4GPA/8lIp/ObPNDme3cLiJL8a8OHsXPSrsJfx7gZzO9oZoY32lPzTRmWVKNMcZkWfWRMcaYLAsKxhhjsiwoGGOMybKgYIwxJsuCgjHGmCwLCsYYY7IsKBhjjMmyoGCMMSbr/wdI180kyfbZxgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.regplot(x = \"Glucose\", y = \"Outcome\", data = diabetes, logistic = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's assess this model by computing the confusion matrix:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[443., 57.],\n", " [138., 130.]])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "confusion_matrix = logit_model.pred_table()\n", "confusion_matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What type of errors are most likely?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the sensitivity and specificity of our model:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sensitivity: 0.886\n", "Specificity: 0.48507462686567165\n" ] } ], "source": [ "true_pos = confusion_matrix[0][0]\n", "false_pos = confusion_matrix[1][0]\n", "false_neg = confusion_matrix[0][1]\n", "true_neg = confusion_matrix[1][1]\n", "sensitivity = true_pos/(true_pos + false_neg)\n", "specificity = true_neg/(true_neg + false_pos)\n", "\n", "print(\"Sensitivity:\",sensitivity)\n", "print(\"Specificity:\",specificity)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Based on the sensitivity and specificity, what are the strong and weak points of the model?\n", "\n", "Remember that by default the confusion matrix uses 0.5 as the cut-off for whether a y value indicates 0 or 1. We can change that by passing in the new cut-off as a parameter. For example, to interpret values >= 0.7 as 1 and < 0.7 as 0, use the code:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[484., 16.],\n", " [195., 73.]])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "confusion_matrix = logit_model.pred_table(0.7)\n", "confusion_matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How did the confusion matrix change? Do you think this new model is better or worse? Recompute the sensitivity and specificity." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sensitivity: 0.968\n", "Specificity: 0.27238805970149255\n" ] } ], "source": [ "true_pos = confusion_matrix[0][0]\n", "false_pos = confusion_matrix[1][0]\n", "false_neg = confusion_matrix[0][1]\n", "true_neg = confusion_matrix[1][1]\n", "sensitivity = true_pos/(true_pos + false_neg)\n", "specificity = true_neg/(true_neg + false_pos)\n", "\n", "print(\"Sensitivity:\",sensitivity)\n", "print(\"Specificity:\",specificity)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How did the sensitivity and specificity change? Are they better or worse?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Logistic Regression with Multiple Independent Variables\n", "\n", "Let's compute a logistic regression model using all of the columns as independent variables." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do any of the variables have higher p-values? If so, let's create a new logistic regression model without those columns." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What's the equation for this logistic regression model?\n", "\n", "As before, let's compute the confusion matrix:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do you think this is an improvement over the model based only on glucose?\n", "\n", "Compute the sensitivity and specificity:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How do the sensitivity and specificity compare to the simpler model? Play around with the cut-off to see what number gives the best results.\n", "\n", "### Challenges\n", "- To better understand the data, plot the distributions of the Glucose column for people with diabetes and without diabetes as overlapping histograms. How does this graph compare to the scatterplot of glucose vs. outcome? Which gives more information?\n", "- (Very challenging) A *Receiver Operating Characteristic curve* or *ROC curve* gives information about the trade-off between the true positive rate (sensitivity) and the false positive rate (1 - specificity) as the cut-off changes. To plot such a curve, use a loop to compute the true positive rate and false positive rate for multiple cut-offs (eg. 0.1, 0.2, ..., 0.8, 0.9), and plot a line plot of these values with the false positive rate on the x axis and the true positive rate on the y axis." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.4.8" } }, "nbformat": 4, "nbformat_minor": 2 }